Script to analyze larval size, symbiont density, and examine correlations between physiological responses.
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# Larval size data
size <- read_csv("Mcap2020/Data/Physiology/Size/larval_size.csv")
#metadata
metadata <- read_csv("Mcap2020/Data/lifestage_metadata.csv")
size <- left_join(size, metadata)
size$hpf <- as.factor(size$hpf)
Prep data frame.
# Calculate mean counts for each sample
size <- size %>%
dplyr::select(tube.ID, lifestage, replicate, `area (mm)`, hpf, group)%>%
drop_na()%>% #remove na's that could not be measured
rename(area=`area (mm)`) #rename column
size$lifestage<-as.factor(size$lifestage)
Plot data with mean and standard error for each lifestage.
size %>%
ggplot(aes(x = lifestage, y = area, color = lifestage)) +
labs(x = "",y = "Mean Larval Size (mm^2)") +
geom_jitter(width = 0.1) + # Plot all points
stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1), # Plot standard error
geom = "errorbar", color = "black", width = 0.5) +
stat_summary(fun = mean, geom = "point", color = "black") + # Plot mean
theme_classic()
Present means and standard error of each group and save summary table
size%>%
group_by(lifestage, hpf)%>%
summarise(n=length(area),
Mean=format(round(mean(area), 3), 3),
SE=format(round(sd(area)/sqrt(length(area)),3),3))%>%
rename(Lifestage=lifestage, HPF=hpf)%>%
kbl(caption="Descriptive statistics of larval size across ontogeny")%>%
kable_classic(full_width=FALSE, html_font="Arial")%>%
row_spec(0, bold = TRUE)
| Lifestage | HPF | n | Mean | SE |
|---|---|---|---|---|
| Egg Fertilized | 1 | 40 | 0.19 | 0.004 |
| Embryo 1 | 5 | 40 | 0.184 | 0.003 |
| Larvae 1 | 38 | 80 | 0.156 | 0.007 |
| Larvae 2 | 65 | 40 | 0.158 | 0.003 |
| Larvae 3 | 93 | 40 | 0.184 | 0.003 |
| Larvae 4 | 163 | 40 | 0.193 | 0.004 |
| Larvae 5 | 183 | 40 | 0.196 | 0.005 |
| Larvae 6 | 231 | 40 | 0.312 | 0.014 |
| Recruit 1 | 183 | 40 | 0.193 | 0.006 |
| Recruit 2 | 231 | 33 | 0.27 | 0.021 |
#need to output to csv
size%>%
group_by(lifestage, hpf)%>%
summarise(n=length(area),
Mean=format(round(mean(area), 3), 3),
SE=format(round(sd(area)/sqrt(length(area)),3),3))%>%
rename(Lifestage=lifestage, HPF=hpf)%>%
write_csv(., "Mcap2020/Output/Physiology/larval_size_table.csv")
Plot data as a scatterplot
size$hpf<-as.factor(size$hpf)
size_plot<-size %>%
ggplot(., aes(x = hpf, y = area)) +
#geom_boxplot(outlier.size = 0) +
geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(0.1)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
ylab(expression(bold(paste("Planar Size (mm"^2, ")"))))+
ylim(0,1)+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); size_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
size_plot2<-size %>%
ggplot(., aes(x = hpf, y = area)) +
geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
geom_point(aes(fill=group), pch = 21, size=2, position = position_jitterdodge(0.1)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
ylab(expression(bold(paste("Planar Size (mm"^2, ")"))))+
ylim(0, 0.8)+
theme_classic() +
#geom_text(label="A", x=1, y=2500, size=4, color="black")+ #egg
#geom_text(label="A", x=2, y=2500, size=4, color="black")+ #embryo 1
#geom_text(label="A", x=3, y=2500, size=4, color="black")+ #larvae 1
#geom_text(label="AB", x=4, y=4100, size=4, color="black")+ #larvae 2
#geom_text(label="AB", x=5, y=4100, size=4, color="black")+ #larvae 3
#geom_text(label="AB", x=6, y=4100, size=4, color="black")+ #larvae 4
#geom_text(label="BC", x=6.8, y=4500, size=4, color="black")+ #larvae 5
#geom_text(label="CD", x=7.2, y=6500, size=4, color="black")+ #recruit 1
#geom_text(label="D", x=7.8, y=6500, size=4, color="black")+ #larvae6
#geom_text(label="D", x=8.2, y=8700, size=4, color="black")+ #recruit2
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14)
); size_plot2
Run lmer on cells per larvae by sampling point, specified by sequence of samples taken (life stage, hpf). Use tube ID as random effect.
size_model<-lmer(area~lifestage + (1|tube.ID), data=size)
summary(size_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: area ~ lifestage + (1 | tube.ID)
## Data: size
##
## REML criterion at convergence: -1223.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7656 -0.4461 -0.0426 0.3090 5.9137
##
## Random effects:
## Groups Name Variance Std.Dev.
## tube.ID (Intercept) 0.000000 0.00000
## Residual 0.002973 0.05453
## Number of obs: 433, groups: tube.ID, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.190150 0.008621 423.000000 22.056 < 2e-16 ***
## lifestageEmbryo 1 -0.005650 0.012192 423.000000 -0.463 0.64332
## lifestageLarvae 1 -0.034300 0.010559 423.000000 -3.248 0.00125 **
## lifestageLarvae 2 -0.031825 0.012192 423.000000 -2.610 0.00937 **
## lifestageLarvae 3 -0.006350 0.012192 423.000000 -0.521 0.60277
## lifestageLarvae 4 0.003025 0.012192 423.000000 0.248 0.80417
## lifestageLarvae 5 0.006000 0.012192 423.000000 0.492 0.62290
## lifestageLarvae 6 0.122175 0.012192 423.000000 10.021 < 2e-16 ***
## lifestageRecruit 1 0.003225 0.012192 423.000000 0.265 0.79152
## lifestageRecruit 2 0.079395 0.012823 423.000000 6.192 1.41e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lfstE1 lfstL1 lfstL2 lfstL3 lfstL4 lfstL5 lfstL6 lfstR1
## lfstgEmbry1 -0.707
## lifestgLrv1 -0.816 0.577
## lifestgLrv2 -0.707 0.500 0.577
## lifestgLrv3 -0.707 0.500 0.577 0.500
## lifestgLrv4 -0.707 0.500 0.577 0.500 0.500
## lifestgLrv5 -0.707 0.500 0.577 0.500 0.500 0.500
## lifestgLrv6 -0.707 0.500 0.577 0.500 0.500 0.500 0.500
## lifstgRcrt1 -0.707 0.500 0.577 0.500 0.500 0.500 0.500 0.500
## lifstgRcrt2 -0.672 0.475 0.549 0.475 0.475 0.475 0.475 0.475 0.475
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
qqPlot(residuals(size_model))
## [1] 401 411
leveneTest(residuals(size_model)~lifestage, data=size)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 9 13.482 < 2.2e-16 ***
## 423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(size_model, type=2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## lifestage 0.91654 0.10184 9 423 34.253 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Violation in normality and variance assumptions. Conduct non-parametric test (Kruskal Wallis).
kruskal.test(area~lifestage, data=size)
##
## Kruskal-Wallis rank sum test
##
## data: area by lifestage
## Kruskal-Wallis chi-squared = 156.39, df = 9, p-value < 2.2e-16
Significant difference in size between lifestages.
View posthoc comparisons for differences between lifestages.
emm = emmeans(size_model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display
## lifestage emmean SE df lower.CL upper.CL .group
## Larvae 1 0.156 0.00610 9.4 0.142 0.170 A
## Larvae 2 0.158 0.00862 38.0 0.141 0.176 AB
## Larvae 3 0.184 0.00862 38.0 0.166 0.201 AB
## Embryo 1 0.184 0.00862 38.0 0.167 0.202 AB
## Egg Fertilized 0.190 0.00862 38.0 0.173 0.208 AB
## Larvae 4 0.193 0.00862 38.0 0.176 0.211 B
## Recruit 1 0.193 0.00862 38.0 0.176 0.211 B
## Larvae 5 0.196 0.00862 38.0 0.179 0.214 B
## Recruit 2 0.270 0.00951 50.5 0.250 0.289 C
## Larvae 6 0.312 0.00862 38.0 0.295 0.330 D
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 10 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
pairs(emm)
## contrast estimate SE df t.ratio p.value
## Egg Fertilized - Embryo 1 0.00565 0.0122 38.0 0.463 1.0000
## Egg Fertilized - Larvae 1 0.03430 0.0106 21.4 3.248 0.0859
## Egg Fertilized - Larvae 2 0.03182 0.0122 38.0 2.610 0.2480
## Egg Fertilized - Larvae 3 0.00635 0.0122 38.0 0.521 0.9999
## Egg Fertilized - Larvae 4 -0.00302 0.0122 38.0 -0.248 1.0000
## Egg Fertilized - Larvae 5 -0.00600 0.0122 38.0 -0.492 1.0000
## Egg Fertilized - Larvae 6 -0.12218 0.0122 38.0 -10.021 <.0001
## Egg Fertilized - Recruit 1 -0.00323 0.0122 38.0 -0.265 1.0000
## Egg Fertilized - Recruit 2 -0.07940 0.0128 44.2 -6.186 <.0001
## Embryo 1 - Larvae 1 0.02865 0.0106 21.4 2.713 0.2293
## Embryo 1 - Larvae 2 0.02618 0.0122 38.0 2.147 0.5082
## Embryo 1 - Larvae 3 0.00070 0.0122 38.0 0.057 1.0000
## Embryo 1 - Larvae 4 -0.00868 0.0122 38.0 -0.712 0.9993
## Embryo 1 - Larvae 5 -0.01165 0.0122 38.0 -0.956 0.9931
## Embryo 1 - Larvae 6 -0.12782 0.0122 38.0 -10.484 <.0001
## Embryo 1 - Recruit 1 -0.00887 0.0122 38.0 -0.728 0.9991
## Embryo 1 - Recruit 2 -0.08505 0.0128 44.2 -6.626 <.0001
## Larvae 1 - Larvae 2 -0.00248 0.0106 21.4 -0.234 1.0000
## Larvae 1 - Larvae 3 -0.02795 0.0106 21.4 -2.647 0.2557
## Larvae 1 - Larvae 4 -0.03732 0.0106 21.4 -3.535 0.0480
## Larvae 1 - Larvae 5 -0.04030 0.0106 21.4 -3.817 0.0264
## Larvae 1 - Larvae 6 -0.15648 0.0106 21.4 -14.819 <.0001
## Larvae 1 - Recruit 1 -0.03753 0.0106 21.4 -3.554 0.0462
## Larvae 1 - Recruit 2 -0.11370 0.0113 26.6 -10.067 <.0001
## Larvae 2 - Larvae 3 -0.02548 0.0122 38.0 -2.089 0.5456
## Larvae 2 - Larvae 4 -0.03485 0.0122 38.0 -2.858 0.1536
## Larvae 2 - Larvae 5 -0.03782 0.0122 38.0 -3.102 0.0909
## Larvae 2 - Larvae 6 -0.15400 0.0122 38.0 -12.631 <.0001
## Larvae 2 - Recruit 1 -0.03505 0.0122 38.0 -2.875 0.1486
## Larvae 2 - Recruit 2 -0.11122 0.0128 44.2 -8.666 <.0001
## Larvae 3 - Larvae 4 -0.00937 0.0122 38.0 -0.769 0.9987
## Larvae 3 - Larvae 5 -0.01235 0.0122 38.0 -1.013 0.9896
## Larvae 3 - Larvae 6 -0.12853 0.0122 38.0 -10.541 <.0001
## Larvae 3 - Recruit 1 -0.00958 0.0122 38.0 -0.785 0.9984
## Larvae 3 - Recruit 2 -0.08575 0.0128 44.2 -6.681 <.0001
## Larvae 4 - Larvae 5 -0.00298 0.0122 38.0 -0.244 1.0000
## Larvae 4 - Larvae 6 -0.11915 0.0122 38.0 -9.772 <.0001
## Larvae 4 - Recruit 1 -0.00020 0.0122 38.0 -0.016 1.0000
## Larvae 4 - Recruit 2 -0.07637 0.0128 44.2 -5.950 <.0001
## Larvae 5 - Larvae 6 -0.11618 0.0122 38.0 -9.528 <.0001
## Larvae 5 - Recruit 1 0.00278 0.0122 38.0 0.228 1.0000
## Larvae 5 - Recruit 2 -0.07340 0.0128 44.2 -5.719 <.0001
## Larvae 6 - Recruit 1 0.11895 0.0122 38.0 9.756 <.0001
## Larvae 6 - Recruit 2 0.04278 0.0128 44.2 3.333 0.0496
## Recruit 1 - Recruit 2 -0.07617 0.0128 44.2 -5.935 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 10 estimates
# Cell count data
sym_counts <- read_csv("Mcap2020/Data/Physiology/CellDensity/symbiont.counts.csv")
sym_counts <- left_join(sym_counts, metadata)
sym_counts$hpf <- as.factor(sym_counts$hpf)
Calculate cells and normalize to either planar size (eggs through metamorphosed recruits) or surface area (attached recruits)
# Calculate mean counts for each sample
df <- sym_counts %>%
dplyr::select(tube.ID, num.squares, matches("count[1-6]")) %>%
gather("rep", "count", -tube.ID, -num.squares) %>%
group_by(tube.ID, num.squares) %>%
summarise(mean_count = mean(count, na.rm = TRUE))
#match in identifying information
df$lifestage<-sym_counts$lifestage[match(df$tube.ID, sym_counts$tube.ID)]
df$total.volume.ul<-sym_counts$total.volume.ul[match(df$tube.ID, sym_counts$tube.ID)]
df$num.individuals<-sym_counts$num.individuals[match(df$tube.ID, sym_counts$tube.ID)]
df$surface.area<-sym_counts$surface.area[match(df$tube.ID, sym_counts$tube.ID)]
df$hpf<-sym_counts$hpf[match(df$tube.ID, sym_counts$tube.ID)]
df$group<-sym_counts$group[match(df$tube.ID, sym_counts$tube.ID)]
df$lifestage<-as.factor(df$lifestage)
df$group<-as.factor(df$group)
# Normalize counts by homogenat volume and surface area
df <- df %>%
mutate(cells.mL = mean_count * 10000 / num.squares,
cells = cells.mL * (total.volume.ul/1000),
cells.ind = cells / num.individuals,
cells.mm = cells / surface.area)
sym_counts<-df
Plot data with mean and standard error for larvae through metamorphosis (these counts have cells/individual). Plot attached recruits separately, these values are in cells per mm2. We will plot cells per unit surface area for all stages in later analyses.
Display cells per individual.
sym_counts %>%
filter(!group=="Attached Recruit")%>%
ggplot(aes(x = lifestage, y = cells.ind, color = lifestage)) +
labs(x = "",y = "Cell Density per larva") +
geom_jitter(width = 0.1) + # Plot all points
stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1), # Plot standard error
geom = "errorbar", color = "black", width = 0.5) +
stat_summary(fun = mean, geom = "point", color = "black") + # Plot mean
theme_classic()
Display cell density per mm2 in attached recruit plugs. Plug 1 = 48 hps, Plug 2 = 72 hps, Plug 3 = 96 hps
sym_counts %>%
filter(group=="Attached Recruit")%>%
ggplot(aes(x = lifestage, y = cells.mm, color = lifestage)) +
labs(x = "",y = "Cell Density per mm2") +
geom_jitter(width = 0.1) + # Plot all points
#stat_summary(fun.data = mean_cl_normal, fun.args = list(mult = 1), # Plot standard error
#geom = "errorbar", color = "black", width = 0.5) +
#stat_summary(fun.y = mean, geom = "point", color = "black") + # Plot mean
theme_classic()
Present means and standard error of each group and save summary table.
sym_counts%>%
group_by(group, hpf, lifestage)%>%
summarise(n=length(cells.ind),
Mean=format(round(mean(cells.ind), 0), 0),
SE=format(round(sd(cells.ind)/sqrt(length(cells.ind)),0),0))%>%
rename(Lifestage=group, HPF=hpf)%>%
kbl(caption="Descriptive statistics of Symbiodiniaceae cell densities per larva across ontogeny")%>%
kable_classic(full_width=FALSE, html_font="Arial")%>%
row_spec(0, bold = TRUE)
| Lifestage | HPF | lifestage | n | Mean | SE |
|---|---|---|---|---|---|
| Attached Recruit | 183 | Plug 1 | 3 | NA | NA |
| Attached Recruit | 231 | Plug 2 | 2 | NA | NA |
| Attached Recruit | 255 | Plug 3 | 3 | NA | NA |
| Egg | 1 | Egg Fertilized | 4 | 1472 | 125 |
| Embryo | 5 | Embryo 1 | 4 | 1831 | 124 |
| Embryo | 38 | Larvae 1 | 4 | 1371 | 242 |
| Embryo | 65 | Larvae 2 | 4 | 2646 | 238 |
| Larvae | 93 | Larvae 3 | 4 | 2692 | 347 |
| Larvae | 163 | Larvae 4 | 4 | 2848 | 206 |
| Larvae | 183 | Larvae 5 | 4 | 3474 | 134 |
| Larvae | 231 | Larvae 6 | 4 | 5142 | 386 |
| Metamorphosed Recruit | 183 | Recruit 1 | 4 | 4829 | 480 |
| Metamorphosed Recruit | 231 | Recruit 2 | 4 | 6424 | 659 |
sym_counts%>%
group_by(group, hpf, lifestage)%>%
summarise(n=length(cells.mm),
Mean=format(round(mean(cells.mm), 0), 0),
SE=format(round(sd(cells.mm)/sqrt(length(cells.mm)),0),0))%>%
rename(Lifestage=group, HPF=hpf)%>%
kbl(caption="Descriptive statistics of Symbiodiniaceae cell densities per mm2 across ontogeny")%>%
kable_classic(full_width=FALSE, html_font="Arial")%>%
row_spec(0, bold = TRUE)
| Lifestage | HPF | lifestage | n | Mean | SE |
|---|---|---|---|---|---|
| Attached Recruit | 183 | Plug 1 | 3 | 21503 | 3084 |
| Attached Recruit | 231 | Plug 2 | 2 | 30977 | 185 |
| Attached Recruit | 255 | Plug 3 | 3 | 42014 | 3847 |
| Egg | 1 | Egg Fertilized | 4 | NA | NA |
| Embryo | 5 | Embryo 1 | 4 | NA | NA |
| Embryo | 38 | Larvae 1 | 4 | NA | NA |
| Embryo | 65 | Larvae 2 | 4 | NA | NA |
| Larvae | 93 | Larvae 3 | 4 | NA | NA |
| Larvae | 163 | Larvae 4 | 4 | NA | NA |
| Larvae | 183 | Larvae 5 | 4 | NA | NA |
| Larvae | 231 | Larvae 6 | 4 | NA | NA |
| Metamorphosed Recruit | 183 | Recruit 1 | 4 | NA | NA |
| Metamorphosed Recruit | 231 | Recruit 2 | 4 | NA | NA |
#need to output to csv
sym_counts%>%
group_by(group, hpf, lifestage)%>%
summarise(n=length(cells.ind),
Mean=format(round(mean(cells.ind), 0), 0),
SE=format(round(sd(cells.ind)/sqrt(length(cells.ind)),0),0))%>%
rename(group=group, HPF=hpf)%>%
write_csv(., "Mcap2020/Output/Physiology/cell_density_table.csv")
Plot data as a scatterplot
sym_counts$hpf<-as.numeric(as.character(sym_counts$hpf))
symb_plot<-sym_counts %>%
filter(!group=="Attached Recruit")%>%
droplevels()%>%
ggplot(., aes(x = hpf, y = cells.ind)) +
#geom_boxplot(outlier.size = 0) +
geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(5)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
ylim(0,9000)+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); symb_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ATTACHED: #BA55D3
Plot data as box plot
symb_plot2<-sym_counts %>%
filter(!group=="Attached Recruit")%>%
droplevels()%>%
ggplot(., aes(x = as.factor(hpf), y = cells.ind)) +
geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
#geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=group), pch = 21, size=4, position = position_jitterdodge(0.2)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"), guide="none")+
ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
ylim(0,9000)+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14)
); symb_plot2
Run ANOVA on cells per larvae by sampling point, specified by sequence of samples taken (life stage, hpf).
sym_ind_model_data<-sym_counts%>%
filter(!group=="Attached Recruit")%>%
droplevels()
sym_ind_model<-aov(cells.ind~lifestage, data=sym_ind_model_data)
summary(sym_ind_model)
## Df Sum Sq Mean Sq F value Pr(>F)
## lifestage 9 102932981 11436998 25.06 1.34e-11 ***
## Residuals 30 13689458 456315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(residuals(sym_ind_model))
## [1] 33 30
leveneTest(residuals(sym_ind_model)~lifestage, data=sym_ind_model_data)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 9 1.7683 0.1166
## 30
Both normality and homogeneity of variance pass.
There is a significant effect of lifestage on cell densities. View posthoc comparisons for differences between lifestages.
emm = emmeans(sym_ind_model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display
## lifestage emmean SE df lower.CL upper.CL .group
## Larvae 1 1371 338 30 681 2061 A
## Egg Fertilized 1472 338 30 782 2161 A
## Embryo 1 1831 338 30 1141 2521 A
## Larvae 2 2646 338 30 1956 3335 AB
## Larvae 3 2692 338 30 2002 3382 AB
## Larvae 4 2848 338 30 2158 3538 AB
## Larvae 5 3474 338 30 2784 4163 BC
## Recruit 1 4829 338 30 4139 5519 CD
## Larvae 6 5142 338 30 4452 5831 D
## Recruit 2 6424 338 30 5734 7114 D
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 10 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
pairs(emm)
## contrast estimate SE df t.ratio p.value
## Egg Fertilized - Embryo 1 -359.6 478 30 -0.753 0.9988
## Egg Fertilized - Larvae 1 100.4 478 30 0.210 1.0000
## Egg Fertilized - Larvae 2 -1174.0 478 30 -2.458 0.3298
## Egg Fertilized - Larvae 3 -1220.4 478 30 -2.555 0.2816
## Egg Fertilized - Larvae 4 -1376.3 478 30 -2.881 0.1554
## Egg Fertilized - Larvae 5 -2002.0 478 30 -4.191 0.0073
## Egg Fertilized - Larvae 6 -3669.9 478 30 -7.683 <.0001
## Egg Fertilized - Recruit 1 -3357.6 478 30 -7.029 <.0001
## Egg Fertilized - Recruit 2 -4952.2 478 30 -10.368 <.0001
## Embryo 1 - Larvae 1 460.0 478 30 0.963 0.9923
## Embryo 1 - Larvae 2 -814.4 478 30 -1.705 0.7842
## Embryo 1 - Larvae 3 -860.8 478 30 -1.802 0.7288
## Embryo 1 - Larvae 4 -1016.7 478 30 -2.128 0.5232
## Embryo 1 - Larvae 5 -1642.4 478 30 -3.438 0.0469
## Embryo 1 - Larvae 6 -3310.3 478 30 -6.930 <.0001
## Embryo 1 - Recruit 1 -2998.0 478 30 -6.276 <.0001
## Embryo 1 - Recruit 2 -4592.6 478 30 -9.615 <.0001
## Larvae 1 - Larvae 2 -1274.4 478 30 -2.668 0.2316
## Larvae 1 - Larvae 3 -1320.8 478 30 -2.765 0.1940
## Larvae 1 - Larvae 4 -1476.7 478 30 -3.091 0.1013
## Larvae 1 - Larvae 5 -2102.4 478 30 -4.401 0.0042
## Larvae 1 - Larvae 6 -3770.3 478 30 -7.893 <.0001
## Larvae 1 - Recruit 1 -3458.0 478 30 -7.239 <.0001
## Larvae 1 - Recruit 2 -5052.6 478 30 -10.578 <.0001
## Larvae 2 - Larvae 3 -46.4 478 30 -0.097 1.0000
## Larvae 2 - Larvae 4 -202.3 478 30 -0.423 1.0000
## Larvae 2 - Larvae 5 -828.0 478 30 -1.733 0.7685
## Larvae 2 - Larvae 6 -2495.9 478 30 -5.225 0.0005
## Larvae 2 - Recruit 1 -2183.6 478 30 -4.571 0.0027
## Larvae 2 - Recruit 2 -3778.2 478 30 -7.910 <.0001
## Larvae 3 - Larvae 4 -155.9 478 30 -0.326 1.0000
## Larvae 3 - Larvae 5 -781.6 478 30 -1.636 0.8201
## Larvae 3 - Larvae 6 -2449.5 478 30 -5.128 0.0006
## Larvae 3 - Recruit 1 -2137.2 478 30 -4.474 0.0035
## Larvae 3 - Recruit 2 -3731.8 478 30 -7.813 <.0001
## Larvae 4 - Larvae 5 -625.8 478 30 -1.310 0.9435
## Larvae 4 - Larvae 6 -2293.7 478 30 -4.802 0.0014
## Larvae 4 - Recruit 1 -1981.3 478 30 -4.148 0.0082
## Larvae 4 - Recruit 2 -3575.9 478 30 -7.486 <.0001
## Larvae 5 - Larvae 6 -1667.9 478 30 -3.492 0.0415
## Larvae 5 - Recruit 1 -1355.6 478 30 -2.838 0.1690
## Larvae 5 - Recruit 2 -2950.2 478 30 -6.176 <.0001
## Larvae 6 - Recruit 1 312.3 478 30 0.654 0.9996
## Larvae 6 - Recruit 2 -1282.3 478 30 -2.684 0.2249
## Recruit 1 - Recruit 2 -1594.6 478 30 -3.338 0.0590
##
## P value adjustment: tukey method for comparing a family of 10 estimates
Output data to file.
sym_counts %>%
write_csv(., file = "Mcap2020/Output/Physiology/calculated_densities.csv")
First, test for correlation between symbiont cell density and larval size to see if there is a relationship.
Generate data frame with summarised size and cell density information for each time point from eggs to metamorphosed recruits, because we have data for size and counts for each sample. We do not include attached recruits here yet, becuase we cannot calculate densities per individual.
#read in data frame generated in previous chunk
sym_counts<-sym_counts%>%
dplyr::select(tube.ID, lifestage, group, hpf, cells.ind, cells.mm)
#grab size data
area<-size%>%
group_by(tube.ID)%>%
summarise(mean_area=mean(area, na.rm=TRUE))
area$tube.ID<-as.factor(area$tube.ID)
sym_counts$hpf<-as.factor(sym_counts$hpf)
corr<-left_join(sym_counts, area)
Generate number of symbiont cells per mm^2 area for each tube.
corr<-corr%>%
mutate(counts_area=cells.ind/mean_area)%>%
mutate(counts_area=ifelse(is.na(counts_area), cells.mm, counts_area)) #add attached recruit data already calculated as cells per mm2
Plot correlation between cell counts (cells per individual) and size (area mm^2).
correlation<-corr %>%
filter(!group=="Attached Recruit")%>%
ggplot(., aes(x = mean_area, y = cells.ind)) +
#geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black", fill="gray") +
geom_point(aes(fill=group), pch = 21, size=4) +
xlab(expression(bold(paste("Larval Size (mm"^2, ")")))) +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30"))+
xlab(expression(bold(paste("Individual Size (mm"^2, ")"))))+
ylab(expression(bold(paste("Symbiont cells individual"^-1))))+
#ylim(0, 9000)+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14)
); correlation
Test relationship with a spearman correlation.
cor.test(corr$mean_area, corr$cells.ind, method=c("spearman"))
##
## Spearman's rank correlation rho
##
## data: corr$mean_area and corr$cells.ind
## S = 3704, p-value = 8.651e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6525328
Significant correlation between size and cell counts. r=0.37, p=0.017
Plot cells per mm^2 as a boxplot.
#order for creating a legend for all plots
corr$group <- factor(corr$group, levels = c("Egg", "Embryo", "Larvae", "Metamorphosed Recruit", "Attached Recruit"))
cells_size_plot<-corr %>%
ggplot(., aes(x = hpf, y = counts_area)) +
geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
geom_point(aes(fill=group), pch = 21, size=4, position = position_jitterdodge(0.4)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"), guide="none")+
scale_color_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"))+
ylab(expression(bold(paste("Symbiont cells mm"^-2))))+
#ylim(2000, 35000)+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14)
); cells_size_plot
Plot as linear relationship.
cells_size_plot2<-corr %>%
ggplot(., aes(x = as.numeric(as.character(hpf)), y = counts_area)) +
#geom_boxplot(aes(color=group), outlier.size = 0, lwd=1) +
geom_point(aes(fill=group, group=group), pch = 21, size=4, position = position_jitterdodge(0.4)) +
geom_smooth(method="lm", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"), guide="none")+
scale_color_manual(name="Lifestage", values=c( "#8C510A", "#DFC27D","#80CDC1", "#003C30", "#BA55D3"))+
ylab(expression(bold(paste("Symbiont cells mm"^-2))))+
#ylim(2000, 35000)+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=14),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14)
); cells_size_plot2
Analyze differences in normalized cell counts by timepoint.
model<-corr%>%
#filter(!group=="Attached Recruit")%>%
#droplevels()%>%
aov(counts_area~lifestage, data=.)
qqPlot(residuals(model))
## [1] 41 30
corr%>%
#filter(!group=="Attached Recruit")%>%
#droplevels()%>%
leveneTest(residuals(model)~lifestage, data=.)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 12 1.1717 0.3399
## 35
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## lifestage 12 3.543e+09 295228137 24.02 2.22e-13 ***
## Residuals 35 4.303e+08 12293680
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
View posthoc differences.
emm = emmeans(model, ~ lifestage)
cld(emm, Letters=c(LETTERS)) #letter display
## lifestage emmean SE df lower.CL upper.CL .group
## Egg Fertilized 7746 1753 35 4187 11305 A
## Larvae 1 8720 1753 35 5161 12279 AB
## Embryo 1 9923 1753 35 6364 13483 ABC
## Larvae 4 14732 1753 35 11173 18291 ABCD
## Larvae 3 14756 1753 35 11197 18315 ABCD
## Larvae 6 16510 1753 35 12951 20069 ABCDE
## Larvae 2 16867 1753 35 13308 20426 BCDE
## Larvae 5 17732 1753 35 14173 21291 CDE
## Plug 1 21503 2024 35 17393 25613 DEF
## Recruit 2 23301 1753 35 19742 26860 DEF
## Recruit 1 25137 1753 35 21578 28696 EF
## Plug 2 30977 2479 35 25944 36010 FG
## Plug 3 42014 2024 35 37905 46124 G
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 13 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
pairs(emm)
## contrast estimate SE df t.ratio p.value
## Egg Fertilized - Embryo 1 -2177.8 2479 35 -0.878 0.9994
## Egg Fertilized - Larvae 1 -974.3 2479 35 -0.393 1.0000
## Egg Fertilized - Larvae 2 -9120.9 2479 35 -3.679 0.0365
## Egg Fertilized - Larvae 3 -7010.7 2479 35 -2.828 0.2348
## Egg Fertilized - Larvae 4 -6986.0 2479 35 -2.818 0.2391
## Egg Fertilized - Larvae 5 -9986.6 2479 35 -4.028 0.0149
## Egg Fertilized - Larvae 6 -8764.6 2479 35 -3.535 0.0517
## Egg Fertilized - Plug 1 -13757.4 2678 35 -5.137 0.0007
## Egg Fertilized - Plug 2 -23231.6 3036 35 -7.651 <.0001
## Egg Fertilized - Plug 3 -34268.5 2678 35 -12.797 <.0001
## Egg Fertilized - Recruit 1 -17390.9 2479 35 -7.014 <.0001
## Egg Fertilized - Recruit 2 -15555.2 2479 35 -6.274 <.0001
## Embryo 1 - Larvae 1 1203.5 2479 35 0.485 1.0000
## Embryo 1 - Larvae 2 -6943.1 2479 35 -2.800 0.2467
## Embryo 1 - Larvae 3 -4832.9 2479 35 -1.949 0.7562
## Embryo 1 - Larvae 4 -4808.2 2479 35 -1.939 0.7620
## Embryo 1 - Larvae 5 -7808.8 2479 35 -3.150 0.1239
## Embryo 1 - Larvae 6 -6586.8 2479 35 -2.657 0.3166
## Embryo 1 - Plug 1 -11579.6 2678 35 -4.324 0.0067
## Embryo 1 - Plug 2 -21053.8 3036 35 -6.934 <.0001
## Embryo 1 - Plug 3 -32090.7 2678 35 -11.983 <.0001
## Embryo 1 - Recruit 1 -15213.1 2479 35 -6.136 <.0001
## Embryo 1 - Recruit 2 -13377.4 2479 35 -5.396 0.0003
## Larvae 1 - Larvae 2 -8146.6 2479 35 -3.286 0.0921
## Larvae 1 - Larvae 3 -6036.4 2479 35 -2.435 0.4450
## Larvae 1 - Larvae 4 -6011.7 2479 35 -2.425 0.4512
## Larvae 1 - Larvae 5 -9012.3 2479 35 -3.635 0.0406
## Larvae 1 - Larvae 6 -7790.3 2479 35 -3.142 0.1259
## Larvae 1 - Plug 1 -12783.0 2678 35 -4.773 0.0019
## Larvae 1 - Plug 2 -22257.3 3036 35 -7.330 <.0001
## Larvae 1 - Plug 3 -33294.2 2678 35 -12.433 <.0001
## Larvae 1 - Recruit 1 -16416.6 2479 35 -6.622 <.0001
## Larvae 1 - Recruit 2 -14580.9 2479 35 -5.881 0.0001
## Larvae 2 - Larvae 3 2110.2 2479 35 0.851 0.9996
## Larvae 2 - Larvae 4 2134.9 2479 35 0.861 0.9995
## Larvae 2 - Larvae 5 -865.7 2479 35 -0.349 1.0000
## Larvae 2 - Larvae 6 356.3 2479 35 0.144 1.0000
## Larvae 2 - Plug 1 -4636.5 2678 35 -1.731 0.8683
## Larvae 2 - Plug 2 -14110.7 3036 35 -4.647 0.0027
## Larvae 2 - Plug 3 -25147.6 2678 35 -9.391 <.0001
## Larvae 2 - Recruit 1 -8270.0 2479 35 -3.336 0.0824
## Larvae 2 - Recruit 2 -6434.3 2479 35 -2.595 0.3499
## Larvae 3 - Larvae 4 24.7 2479 35 0.010 1.0000
## Larvae 3 - Larvae 5 -2975.9 2479 35 -1.200 0.9902
## Larvae 3 - Larvae 6 -1753.9 2479 35 -0.707 0.9999
## Larvae 3 - Plug 1 -6746.7 2678 35 -2.519 0.3935
## Larvae 3 - Plug 2 -16220.9 3036 35 -5.342 0.0004
## Larvae 3 - Plug 3 -27257.8 2678 35 -10.179 <.0001
## Larvae 3 - Recruit 1 -10380.2 2479 35 -4.187 0.0098
## Larvae 3 - Recruit 2 -8544.5 2479 35 -3.446 0.0638
## Larvae 4 - Larvae 5 -3000.6 2479 35 -1.210 0.9896
## Larvae 4 - Larvae 6 -1778.6 2479 35 -0.717 0.9999
## Larvae 4 - Plug 1 -6771.3 2678 35 -2.529 0.3881
## Larvae 4 - Plug 2 -16245.5 3036 35 -5.350 0.0004
## Larvae 4 - Plug 3 -27282.5 2678 35 -10.188 <.0001
## Larvae 4 - Recruit 1 -10404.9 2479 35 -4.197 0.0095
## Larvae 4 - Recruit 2 -8569.2 2479 35 -3.456 0.0624
## Larvae 5 - Larvae 6 1222.0 2479 35 0.493 1.0000
## Larvae 5 - Plug 1 -3770.8 2678 35 -1.408 0.9658
## Larvae 5 - Plug 2 -13245.0 3036 35 -4.362 0.0061
## Larvae 5 - Plug 3 -24281.9 2678 35 -9.067 <.0001
## Larvae 5 - Recruit 1 -7404.3 2479 35 -2.986 0.1733
## Larvae 5 - Recruit 2 -5568.6 2479 35 -2.246 0.5671
## Larvae 6 - Plug 1 -4992.8 2678 35 -1.864 0.8038
## Larvae 6 - Plug 2 -14467.0 3036 35 -4.764 0.0020
## Larvae 6 - Plug 3 -25503.9 2678 35 -9.524 <.0001
## Larvae 6 - Recruit 1 -8626.3 2479 35 -3.479 0.0591
## Larvae 6 - Recruit 2 -6790.6 2479 35 -2.739 0.2753
## Plug 1 - Plug 2 -9474.2 3201 35 -2.960 0.1826
## Plug 1 - Plug 3 -20511.1 2863 35 -7.165 <.0001
## Plug 1 - Recruit 1 -3633.5 2678 35 -1.357 0.9741
## Plug 1 - Recruit 2 -1797.8 2678 35 -0.671 1.0000
## Plug 2 - Plug 3 -11036.9 3201 35 -3.448 0.0636
## Plug 2 - Recruit 1 5840.7 3036 35 1.924 0.7712
## Plug 2 - Recruit 2 7676.4 3036 35 2.528 0.3884
## Plug 3 - Recruit 1 16877.6 2678 35 6.302 <.0001
## Plug 3 - Recruit 2 18713.3 2678 35 6.988 <.0001
## Recruit 1 - Recruit 2 1835.7 2479 35 0.740 0.9999
##
## P value adjustment: tukey method for comparing a family of 13 estimates
Generate summary table of descriptive statistics.
#need to output to csv
corr%>%
group_by(group, hpf, lifestage)%>%
summarise(n=length(counts_area),
Mean_sym_mm2=format(round(mean(counts_area), 0), 0),
SE=format(round(sd(counts_area)/sqrt(length(counts_area)),0),0))%>%
rename(Lifestage=group, HPF=hpf)%>%
write_csv(., "Mcap2020/Output/Physiology/normalized_size_cells_summary.csv")
Generate physiology panel with all variables of interest.
# extract the legend from one of the plots
legend <- get_legend(
# create some space to the left of the legend
cells_size_plot + theme(legend.box.margin = margin(1,1,1,1))
)
#remove legends from plots
size_plot2<-size_plot2+theme(legend.position="none")
symb_plot2<-symb_plot2+theme(legend.position="none")
r_corr_plot<-r_corr_plot+theme(legend.position="none")
cells_size_plot_l<-cells_size_plot+theme(legend.position="none")
#assemble plots
all_plots<-plot_grid(size_plot2, symb_plot2, cells_size_plot_l, r_corr_plot, labels = c("A", "B", "C", "D"), label_size=18, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,0.8), align="h")
all_plots_legend<-plot_grid(all_plots, legend, rel_widths = c(4, 0.5), ncol=2, nrow=1)
ggsave(file="Mcap2020/Figures/Physiology/Physiology_figure.pdf", all_plots_legend, dpi=300, width=24, height=6, units="in")
ggsave(file="Mcap2020/Figures/Physiology/Physiology_figure.png", all_plots_legend, dpi=300, width=24, height=6, units="in")
Early life history physiology